*! by Hung-Jen Wang
*!    Department of Economics, National Taiwan University
*!    wangh@ntu.edu.tw

*! Generate random draws from truncated standard normal distribution.
*! Update on 24July00 to allow for variables in Left or Right
*! Update on 17Nov00 to allow for non-standard normal distribution.

capture program drop gentrun
program define gentrun
   version 6
   syntax newvarlist [if] [in] [, DIST(string) Left(string) Right(string) SCALAR]

   tempvar tem1 tem2 mltrn mrtrn myvar mmean mvari


   if "`left'" != "" & "`right'" != "" {
      if `left' >= `right' {
         di in red /*
    */ "The left truncation point cannot be >= the right truncation point."
         exit 198
      }
   }


   if "`dist'" != "" { /* for a possibly non-standard normal distribution */
      tokenize `dist'
      if "`2'" == ""{
         di in red "If DIST is used, you must specify 2 numerical numbers, first for mean and second for variance."
         exit 198
      }

      quie gen double `mmean' = `1'
      quie gen double `mvari' = `2'
   }
   else {
      quie gen `mmean' = 0
      quie gen `mvari' = 1
   }


     if "`left'" == ""{
         gen `mltrn' = 0 /* if not truncated at left, the implied CDF begins at 0 */
     }
     else {
         gen double `mltrn' = normprob((`left'-`mmean')/sqrt(`mvari'))
     }

     if "`right'" == ""{
         gen `mrtrn' = 1  /* if not truncated at right, the implied CDF is up to 1 */
     }
     else {
         gen double `mrtrn' = normprob((`right'-`mmean')/sqrt(`mvari'))
     }

    tokenize `varlist'

   if "`scalar'" ~= "" {  /* if create a scalar */
    while "`1'" ~= ""{
      scalar `tem1' = (`mrtrn' - `mltrn')*uniform() + `mltrn'
      scalar `tem2' = invnorm(`tem1')
      scalar `1' = `mmean' + `tem2'*sqrt(`mvari')
      mac shift
    }
   }
   else{          /* if create a variable */
    while "`1'" ~= ""{
        quie gen double `tem1' = (`mrtrn'-`mltrn')*uniform() + `mltrn' `if' `in'
        quie gen double `tem2' = invnorm(`tem1') `if' `in' /* convert CDF back to x */
        quie gen `typlist' `1' = `mmean' + `tem2'*sqrt(`mvari')
        drop `tem1' `tem2'
        mac shift
    }
   }

end
exit

* x -> normprob -> cdf; cdf -> invnorm -> x.
* left, right: x
* mltrn, mrtrn: cdf
*       scalar `tem1' = (`mrtrn' - `mltrn')*uniform() + `mltrn'
* The logic of the above line is this: The (mrtrn - mltrn) is the CDF
* after the truncation.  For example, if right is not truncated and left
* is truncated at 0, so I am looking for the right-half distribution,
* then (mrtrn - mltrn) would equal 0.5, which is the CDF of the right-
* half distribution.  Now I know the rv is going to be generated from
* a distribution whose before-normalized CDF is 0.5.  Since I want a
* random draw, so I pick any z whose corresponding CDF would be in the
* sensible range, which in the above case is in the range between 0 and
* 0.5.  This is what the "uniform()" does: randomly pick a CDF value
* in the sensible range, ie. 0 and 0.5.  But note that when I am doing
* this, I am picking a point on the "left-half" of the distribution,
* rather than the right-half.  This is because CDF between 0 and 0.5
* is actually counted from the left.  Since what I really want is
* the right-half, I add in "mltrn", which is the left-truncated CDF, so
* it moves the point to the right and to the appropriate place.
*       scalar `1' = invnorm(`tem1')
* Now that I have the Z's corresponding CDF, I simply use invnorm to
* put it back to the Z value.


* About random draws from non-standard normal.
* The thing is simple.  If x is drew from N(mu, sig2), then let
* z = (x-mu)/sig and draw it from N(0, 1).  After that, convert
* it back by x = mu + sig*z.  Operationally, I have to convert the
* truncation points as well.
